home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / ada / gnat1792.zip / gnat179b / t-adainc / a-ngelfu.adb < prev    next >
Text File  |  1994-05-19  |  20KB  |  763 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT RUNTIME COMPONENTS                          --
  4. --                                                                          --
  5. --                     A D A . N U M E R I C S . G E F                      --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.2 $                              --
  10. --                                                                          --
  11. --           Copyright (c) 1992,1993,1994 NYU, All Rights Reserved          --
  12. --                                                                          --
  13. -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14. -- terms of the  GNU General Public License as published  by the Free Soft- --
  15. -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16. -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17. -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18. -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19. -- for  more details.  You should have  received  a copy of the GNU General --
  20. -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21. -- to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. --
  22. --                                                                          --
  23. ------------------------------------------------------------------------------
  24.  
  25. --  This body is specifically for using an Ada interface to C math.h
  26. --  to get the computation engine
  27. --  Many special cases are handled locally to avoid unnecessary calls
  28. --  All Ada required exception handling is provided.
  29. --  This is not a "strict" implementation.
  30. --  This is an early, trial, implementation; simplified for early gnat.
  31. --  uses: sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
  32. --  sinh, cosh, tanh
  33.  
  34. with Ada.Numerics.Aux; -- interface to C library via math.h is in this package
  35.  
  36.  
  37. package body Ada.Numerics.Generic_Elementary_Functions is
  38.  
  39.    subtype Double is Ada.Numerics.Aux.Double;
  40.  
  41.    Pi      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41972;
  42.    Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  43.    Two_Pi  : Float_Type'Base;
  44.    Half_Pi : Float_Type'Base;
  45.    Fourth_Pi : Float_Type'Base;
  46.    Epsilon : Float_Type'Base;
  47.    Square_Root_Epsilon : Float_Type'Base;
  48.    Half_Log_Epsilon : Float_Type'Base;
  49.    Log_Last : Float_Type'Base;
  50.    Log_Inverse_Epsilon : Float_Type'Base;
  51.  
  52.    function Exact_Remainder
  53.      (X : Float_Type'Base;
  54.       Y : Float_Type'Base)
  55.       return Float_Type'Base
  56.    is
  57.       Denominator : Float_Type'Base := abs X;
  58.       Divisor : Float_Type'Base := abs Y;
  59.       Reducer : Float_Type'Base;
  60.       Sign : Float_Type'Base := 1.0;
  61.    begin
  62.       if Y = 0.0 then
  63.          raise Constraint_Error;
  64.       elsif X = 0.0 then
  65.          return 0.0;
  66.       elsif X = Y then
  67.          return 0.0;
  68.       elsif Denominator < Divisor then
  69.          return X;
  70.       end if;
  71.  
  72.       while Denominator >= Divisor loop
  73.          --  put divisors mantissa with denominators exponent to make reducer
  74.          Reducer := Divisor;
  75.  
  76.          begin
  77.             while Reducer * 1_048_576.0 < Denominator loop
  78.                Reducer := Reducer * 1_048_576.0;
  79.             end loop;
  80.          exception
  81.             when others => null;
  82.          end;
  83.  
  84.          begin
  85.             while Reducer * 1_024.0 < Denominator loop
  86.                Reducer := Reducer * 1_024.0;
  87.             end loop;
  88.          exception
  89.             when others => null;
  90.          end;
  91.  
  92.          begin
  93.             while Reducer * 2.0 < Denominator loop
  94.                Reducer := Reducer * 2.0;
  95.             end loop;
  96.          exception
  97.             when others => null;
  98.          end;
  99.  
  100.          Denominator := Denominator - Reducer;
  101.       end loop;
  102.  
  103.       if X < 0.0 then
  104.          return -Denominator;
  105.       else
  106.          return Denominator;
  107.       end if;
  108.    end Exact_Remainder;
  109.  
  110.    ----------------
  111.    -- Local_Atan --
  112.    ----------------
  113.  
  114.    function Local_Atan
  115.      (Y : Float_Type'Base;
  116.       X : Float_Type'Base := 1.0)
  117.       return Float_Type'Base
  118.    is
  119.       Z : Float_Type'Base;
  120.       Raw_Atan : Float_Type'Base;
  121.    begin
  122.       if abs Y > abs X then
  123.          Z := abs (X / Y);
  124.       else
  125.          Z := abs (Y / X);
  126.       end if;
  127.  
  128.       if Z < Square_Root_Epsilon then
  129.          Raw_Atan := Z;
  130.       elsif Z = 1.0 then
  131.          Raw_Atan := Pi / 4.0;
  132.       elsif Z < Square_Root_Epsilon then
  133.          Raw_Atan := Z;
  134.       else
  135.          Raw_Atan := Float_Type'Base (Ada.Numerics.Aux.Atan (Double (Z)));
  136.       end if;
  137.  
  138.       if abs Y > abs X then
  139.          Raw_Atan := Half_Pi - Raw_Atan;
  140.       end if;
  141.  
  142.       if X > 0.0 then
  143.          if Y > 0.0 then
  144.             return Raw_Atan;
  145.          else                 --  Y < 0.0
  146.             return -Raw_Atan;
  147.          end if;
  148.  
  149.       else                    --  X < 0.0
  150.          if Y > 0.0 then
  151.             return Pi - Raw_Atan;
  152.          else                  --  Y < 0.0
  153.             return -(Pi - Raw_Atan);
  154.          end if;
  155.       end if;
  156.    end Local_Atan;
  157.  
  158.    ----------
  159.    -- Sqrt --
  160.    ----------
  161.  
  162.    function Sqrt (X : Float_Type'Base) return Float_Type'Base is
  163.    begin
  164.       if X < 0.0 then
  165.          raise Argument_Error;
  166.       elsif X = 0.0 then
  167.          return X; -- may be +0.0 or -0.0
  168.       elsif X = 1.0 then
  169.          return 1.0; -- needs to be exact for good COMPLEX accuracy
  170.       end if;
  171.  
  172.       return Float_Type'Base (Ada.Numerics.Aux.Sqrt (Double (X)));
  173.    end Sqrt;
  174.  
  175.    -------------------------
  176.    --  Log (natural base) --
  177.    -------------------------
  178.  
  179.    function Log (X : Float_Type'Base) return Float_Type'Base is
  180.    begin
  181.       if X < 0.0 then
  182.          raise Argument_Error;
  183.       elsif X = 0.0 then
  184.          raise Constraint_Error;
  185.       elsif X = 1.0 then
  186.          return 0.0;
  187.       end if;
  188.  
  189.       return Float_Type'Base (Ada.Numerics.Aux.Log (Double (X)));
  190.    end Log;
  191.  
  192.    --------------------------
  193.    -- Log (arbitrary base) --
  194.    --------------------------
  195.  
  196.    function Log (X, Base : Float_Type'Base) return Float_Type'Base is
  197.    begin
  198.       if X < 0.0 then
  199.          raise Argument_Error;
  200.       elsif Base <= 0.0 or else Base = 1.0 then
  201.          raise Argument_Error;
  202.       elsif X = 0.0 then
  203.          raise Constraint_Error;
  204.       elsif X = 1.0 then
  205.          return 0.0;
  206.       end if;
  207.  
  208.       return Float_Type'Base (Ada.Numerics.Aux.Log (Double (X))
  209.                       / Ada.Numerics.Aux.Log (Double (Base)));
  210.    end Log;
  211.  
  212.    ---------
  213.    -- Exp --
  214.    ---------
  215.  
  216.    function Exp (X : Float_Type'Base) return Float_Type'Base is
  217.       Result : Float_Type'Base;
  218.    begin
  219.       if X = 0.0 then
  220.          return 1.0;
  221.       elsif X > Log_Last then
  222.          raise Constraint_Error;
  223.       end if;
  224.  
  225.       return Float_Type'Base (Ada.Numerics.Aux.Exp (Double (X)));
  226.    exception
  227.       when others =>
  228.          raise Constraint_Error;
  229.    end Exp;
  230.  
  231.    ----------
  232.    -- "**" --
  233.    ----------
  234.  
  235.    function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
  236.       Result : Float_Type'Base;
  237.    begin
  238.       if Left = 0.0
  239.         and then Right = 0.0
  240.       then
  241.          raise Argument_Error;
  242.       elsif Left < 0.0 then
  243.          raise Argument_Error;
  244.       elsif Right = 0.0 then
  245.          return 1.0;
  246.       elsif Left = 0.0 then
  247.  
  248.          if Right < 0.0 then
  249.             raise Constraint_Error;
  250.          else
  251.             return 0.0;
  252.          end if;
  253.  
  254.       elsif Left = 1.0 then
  255.          return 1.0;
  256.       elsif Right = 1.0 then
  257.          return Left;
  258.       elsif Right = 2.0 then
  259.          return Left * Left;
  260.       end if;
  261.  
  262.       return Float_Type'Base
  263.                (Ada.Numerics.Aux.pow (Double (Left), Double (Right)));
  264.    exception
  265.       when others =>
  266.          raise Constraint_Error;
  267.    end "**";
  268.  
  269.  
  270.    -------------------------
  271.    -- Sin (natural cycle) --
  272.    --------------